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Calculations of Raman scattering intensities in spin 1/2 square- lattice Heisenberg model, using 
the Fleury-Loudon-Elliott theory, have so far been unable to describe the broad line shape and 
asymmetry of the two magnon peak found experimentally in the cuprate materials. Even more 
notably, the polarization selection rules are violated with respect to the Fleury-Loudon-Elliott the- 
ory. There is comparable scattering in B\ g and Ai g geometries, whereas the theory would predict 
scattering in only Bi g geometry. We review various suggestions for this discrepency and suggest 
that at least part of the problem can be addressed by modifying the effective Raman Hamiltonian, 
allowing for two-magnon states with arbitrary total momentum. Such an approach based on the 
Sawatzsky-Lorenzana theory of optical absorption assumes an important role of phonons as momen- 
tum sinks. It leaves the low energy physics of the Heisenberg model unchanged but substantially 
alters the Raman line-shape and selection rules, bringing the results closer to experiments. 
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I. INTRODUCTION 

00 ; 

I Raman scattering provides an important tool for examining the structure of antiferromagnetic materials. Even 
(f) , though optical processes in Mott insulators necessarily depend on the energy-bands in a complex way, the Fleury- 
Loudon-Elliott theoryM allows one to bypass that complexity and develop a theory for the line-shape of the Raman 
spectra entirely within the framework of pure spin models. This theory has been highly successful in many cases and 
is the primary reason why the Raman scattering becomes an investigative tool for these class of materials. However, in 
case of the cuprates such as La2Cu04, the stochiometric parent compounds of the high temperature superconducting 
materials, the Fleury-Loudon-Elliott theory runs into several difficulties. This has been a subject of intense debate, 
and many explanations have been proposed ranging from the inadequacy of the theory to novel and exotic microscopic 
. physics in these materials. The goal of this paper is to review the various explanations and to examine how far a simple 
modification to the effective "Raman Hamiltonian" that allows two-magnon states with arbitrary total momentum, 
O ■ can help bridge the gap between theory and experiments. 

O ■ R is fair to say, that Raman scattering provided some of the earliest accurate estimates for the antiferromagnetic 
\ exchange constant in the cupratesD This, by itself, is proof that the peak-frequency of the Raman scattering intensity 
matches reasonably with theoretical expectations. Furthermore- the lineshape of the spectra is reasonably universal 
r> | from one material to another within the insulating cuprates pa and although there are details in the shape whose 
dependence on the incident photon energy can be clearly recognized, the gross features of the lineshape are largely 
independent of such resonance effects. Thus, the discrepency between theory and experiment only come in when a 
more detailed calculation of the lineshape is performed within the Fleury-Loudon-Elliott theory. The experimental 
spectra is much broader, perhaps by about a factor of 3, and has a clear asymmetry extending towards high energies. 
The most notable discrepency with the theory is that in the experiments there is comparable scattering intensity 
in Ai g and B\ g polarizations of incident and outgoing light, whereas theory predicts scattering predominantly in 
B\g geometry only. The fact that the main features of the spectra are so universal suggests that it is intrinsic and 
significant. 

The theoretical work focusing on these discrepencies can be grouped into the following categories: 
(i) Inaccuracies of numerical calculations: Even given a system well described by a nearest-neighbor Heisenberg 
model and an effective spin-Hamiltonian which describes the Raman scattering process, the calculation of Raman 
scattering lineshape remains a challenging task. The spin-wave theory.^ which works well for higher spin and higher 
dimensional systeps, need not be accurate for a 2D spin-half system. Improved calculations have involved higher-order 
spin-wave theory,£9 series-expansions, til, exact-diagonalization of small systemaia and finite temperature Quantum 
Monte Carlo simulations £3 These calculations have established the first few moments of the spectra quite well. The 
Quantum Monte Carlo calculation is perhaps the best in terms of getting the lineshape correct and suggests that the 
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actual lineshape can be fairly different from spin- wave theory. It maybe both broader than spin- wave theory and have 
some of the high energy asymmetry, but perhaps not as much as in the experiments. 

(ii) The Heisenberg model is not good enough for the cuprates: Other work has focused on extending the nearest- 
neighbor Heisenberg model in order to get better agreement with the experiments. Foi^xample, one could introduce 
second neighbor antiferromagnetic interactions to explain scatteringpirL-jdi^ geometries. til A more radical proposal has 
been the possibility of substantial or dominant ring-exchange termsJi-TcJ which can dramatically broaden the spectra. 
The consistency of such an approach with other measurements (.most notably neutron scattering) has not been shown. 

(iii) Lineshape depends on resoance: Chubukov and Frenkeltj ancLindependently Kampf et al. have argued that 
the lineshape does depend on the frequency of incident photon energyL3 and these features can also make the spectra 
appear broader and give enhanced scattering at higher frequencies. 

(iv) Other degrees of freedom, most notably phonons are important: It has been argued by several authors that 
coupling between spin and phonons can lead to substantial broadening of the spectra. Calculations in this respect 
have included modeling phonons by substantial modulation of local coupling constantsEj as well as by spin-wave 
theory.Ea Again, the consistency of strong spin-phonon couplings to neutron scattering and other measurements have 
not been shown. In particular, the fact that neutron scattering measurements, especially the tempeiajbure dependent 
correlation length £(T) and the spin-dynamics, agree remarkably well with the Heisenberg mode£!rE3 does not leave 
room for such couplings. 

(v) Magnons are not good elementary excitations at short wavelengths: One of the most exciting suggestions from 
a physics point of vieiM-hap been to invoke spinons and not magnons as elementary excitations of the system, at least 
at short wavelengths. c5c3 Such an approach naturally leads to a much broader spectra, and can be considered to be 
successful at a phenomenological level. The primary difficulty with this approach is that the existence of spinon-likc 
excitations in two-dimensions remains highly controversial. I— ■ 

(vi) The need to go beyond the spin-subspace to describe the scattering process: The work of Shastry and ShraimanEZl 
has presented a comprehensive theoretical framework for understanding the Fleury-Loudon-Elliott scheme for effective 
Raman Hamiltonians starting from an electronic Hamiltonian. However, the cuprate materials are far from the large-U 
limit where such a scheme can be rigorously shown to work, and thus multiple bands and detailed band-structure may 
play a role here. However, as noted before, the fact that the spectral features are reasonably universal over different 
family of materials suggests a more generic explanation may be appropriate. 

In this paper, we primarily concern ourselves with the numerical calculation of the Raman scattering lineshape 
with a modified effective Hamiltonian. Such an approach does not alter the ground state properties and elementary 
excitations of the system, but only the way in which the Raman scattering process is described within the spin 
subspace. The basic idea is based on the work of Sawatzsky and LorenzanacEl for optical absorption in the cuprates. 
They argued that the optical absorption was assisted by phonons, whose role can be incorporated into the theory by 
simply assuming that they acted as momentum sinks. Thus optical absorption could proceed through excitation of 
two-magnons with arbitrary total momentum. Here we explore the analogous situation for Raman scattering aided by 
phonons. This immedeately leads to scattering in both B lg and A lg geometries. Furthermore, the spectral features 
come closer to experiments. Given our finite-size numerical calculations, it is difficult to say whether these are now 
in complete agreement with experiments. 



II. PHONON ASSISTED SCATTERING: THE SINGLE SITE OPERATOR (SSO) 

Let us first examine the rationale behind the very successful Fleury-Loudon§ theory as formulated by Elliotti. 
Although Raman scattering proceeds through virtual charge excitations, the scattering process can be described by 
an effective spin Hamiltonian, simply by incorporating the important symmetries of the problem. This is possible 
because the initial and final states both lie well below the charge-gap and thus the resulting excitation must be a 
pure spin excitation. Since light has very long-wavelength and the scattering involves the electric field and not the 
magnetic field, the effective Raman Hamiltonian must have zero total momentum and be a spin-singlet. It must be 
linear in the polarizations of incoming and outgoing electric field vectors and must be a scalar. If we further assume 
the dominance of nearest-neighbor superexchange, the effective Raman Hamiltonian, is essentially fully determined 
apart from an overall multiplicative constant. It takes the Fleury-Loudon-Elliott form: 

Hr = ^ fen ■ hj)feout ■ hj)Si ■ Sj, (1) 

<ij> 

where, the sum runs over the nearest-neighbor pairs, Ci n and e ou t are the incoming and outgoing electric field polar- 
ization vectors and fij is a unit vector connecting the sites i and j. 



2 



Thus in the B\ g configuration, where the incoming and outgoing light are polarized in the plane of the copper- 
oxides at right angles to each other and at an angle 45 degrees from the x and y axes of the CuC>2 lattice, the effective 
scattering operator becomes: 

o Blg = ]T Si ( 2 ) 

<ij>,x <ij>,y 

where the first sum is over the nearest neighbor bonds parallel to the x- axis, and the second sum is over the nearest 
neighbor bonds parallel to the y-axis. 

In contrast, the effective Hamiltonian vanishes in the B 2g configuration, where the incoming and outgoing light are 
polarized in the plane of the copper-oxides at right angles to each other, with one being along the x and the other 
along the y axis. In the A\ g configuration, the Fleury-Loudon-Elliott operator becomes 

Al3 = ]T ]T Si -4 (3) 

<ij>,x <ij>-,y 

which is just the Hciscnberg Hamiltonian and, thus, results in no scattering. Thus the theory predicts scattering in 
B\g geometry only. The spectra obtained by treating this effective Hamiltonian in the two-magnon approximation 
for S > 1, provide a remarkably accurate description of the experiments in K2N1F4 and other materials.^ Numerical 
calculations for S = 1/2, and their lack of agreement with the cuprate materials will be discussed in the following 
sections. 

Here we will examine the possibility that the phonons or impurities play an important role in the Raman scattering 
process, even though they do not much effect the system in the .absence of incident light. One way to mimic the 
role of phonons follows from the work of Lorenzana and SawatzkyE3 on optical absorption in antiferromagnets. The 
key effect of phonons, in their theory, is to act as a momentum sink, allowing absorption via two-magnon states 
of arbitrary total momentum™ This theory has proved to be very successful in describing optical absorption in the 
quasi- ID material Sr2CuO3.ES 

We can incorporate this idea of phonons acting as momentum sinks in Raman scattering by modifying the Raman 
scattering operator. The most natural choice is to consider the following single site operator (SSO) for the B\ g 
configuration: 

o Blg = E ( 4 ) 



<j>,x <j>,y 



And, for the Ai g configuration: 



°M 9 = E S? °'^ + E So-Sj. (5) 



<j>,x <j>,y 



Notice that the latter no longer commutes with the antiferromagnetic Heisenberg Hamiltonian, and can thus produce 
scattering in the A\ g channel. In general, one would expect that by including two-magnons at arbitrary momentum the 
two magnons will scatter less with each other in the final state and thus lead to a broadening of the spectra. Whether 
this effect combined with quantum fluctuations can lead to a spectra consistent with the experiments becomes a 
numerical issue. 



III. COMPUTATIONAL METHODS 



Results shown in this paper will be based on exact diagonalization computations of Raman spectra. The antifer- 
romagnetic Heisenberg model is used for systems o£, 16 and 26 sites. To obtain a ground state vector |*0o > having 
energy Eq, a conjugate gradient method was usedHS Once that was accomplished, it was possible to compute zero- 
temperature Raman spectra using a variety of methods, which we will now discuss. Let us assume that our scattering 
operator is O; the equation for the scattering intensity / at the shifted frequency lo has the form 

I(u) = --Im[< MO^ - = ]_ . ^O|^o >], (6) 

it oj + hi® + ie — H 

where H is the Hamiltonian of the system and e is a small real number introduced to allow computation. This 
equation can also be expressed in a Fermi's golden rule form, 
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I(oj) =^2\< i>n\O\i> > \ 2 5{uo - (E n - B )), 

■n 



(7) 



where \ip n > and E n are eigenvectors and eigenvalues of the system. 

There are many possible ways to perform this calculation usingj-these two equation forms. One standard method is 
to use a continued fraction calculation on the first form. Dagottotj describes how this calculation can be performed. 
More recent techniques relying on the second form of the scattering equation are simpler to implement, however. The 
first one we shall examine, which is sometimes called the spectral decoding technique, was first introduced by Loh 
and Campbell.Lil Let us define a set of vectors |0„ > using the well-known Lanczos iteration technique^ 



\(t>i >= H\(f) > |0 O >, (9) 

< <Po |0o > 

and 

U ^ H-U < 0n|g|0„ > , , <0»|0n> ,.,„,. 

n+ l >= > — — \(j) n > — \<Pn-l > ■ (10) 

< <Pn\9n > < <Pn-l\<Pn-l > 

With this set of vectors defined, we now have a simple tridiagonal form for the Hamiltonian matrix that can be easily 
diagonalized. We can now say that the eigenvectors \ijj n > are related to the \<fi n > by the relationship 

l^n >= X>3m >• (11) 

m 

It can now be shown that 

| < Vn|0|^o > I 2 = K\ 2 < ^olOtOlVo > • (12) 

The final spectrum can be displayed by replacing the Dirac delta functions in Eq. (Q) with finite Lorenzians of an 
arbitrary width. 

Spectral decoding is a very useful technique, but it has some disadvantages. It relies on the Lanczos method for 
eigenvector computation above the ground state, and it is known that the Lanczos method can produce eigensolutions 
which are either incorrect or are duplicates of other pSolutions found previously. Techniques exist for checking the 
validity of solutions provided by the Lanczos methodJ23 which we shall call sorting, but they can be cumbersome. It 
would be preferable to use another technique where sorting is not necessary. 

For the spectra computed in this paper, the kernel polynomial method^ (KPM) was used. In KPM, a convergent 
approximation to the true spectrum is computed using Chebyshev polynomials. The delta function in Eq. (0) is 
replaced with a Chebyshev expansion of the delta function, and Gibbs damping factors are included to eliminate the 
Gibbs phenomenon. Calculations are performed using the operator X instead of H, where X is simply rescaled so 
that all energies lie between -1 and 1. Similarly, we use x instead of u>, where x is u> rescaled to lie between and 2. 
The final calculation involved is 

/ ( :E ) = — n \ ,xo =[ffo^o + 2 V g m T m (x - l)/x m ] (13) 

7iyi-(£-ir 

where the T m are Chebyshev polynomials, the g m are Gibbs damping factors, and the moments fi m are defined by 

fi m =< ^o|0 + T ro (X)0|Vo > ■ (14) 

The T m (X) here are Chebyshev polynomials of the operator X . In practice, the moments are most easily calculated 
using Chebyshev recurrence relations. Using these relations, computing M moments requires only 4^ + 1 calculations. 

KPM results are equivalent to those of other methods mentioned above. KPM is used here because it is simpler to 
implement computationally than other methods for a given level of accuracy. 
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IV. RESULTS 



Now let us examine some of our results. Calculations for the 16 site model were performed on a 200MHz personal 
computer. The Nccl state was used as a starting point for ground state calculations. Spin flip symmetry was used 
to reduce the final size of the Hilbert space to 6435. Memory requirements were minimal, and calculations were 
accomplished in minutes. The 26 site system spectra were computed on various machines with Alpha processors. 
Again, a spin flip symmetry was the only symmetry used, the Neel state was used as an initial approximation to 
the ground state, and the Hilbert space had a dimensionality of 5,200,300. Several hundred megabytes of RAM were 
required, and the calculations were completed in several hours time. 

Fig. [I] shows some computed spectra in the B\ g configuration for a 16 site Heisenberg model of the square lattice 
with periodic boundary conditions. If we examine the Fleury-Loudon-Elliott spectrum (the solid curve in part (a)) 
and the SSO spectrum (the solid curve in part (b)) we see that there is much more-activity in the SSO spectrum, 
and that its greatest activity occurs more toward the peak determined by experimentliil for La2Cu04 (dashed curves). 
In Fig. U, we see the SSO spectrum broadened and shifted slightly so that its main peak is in the same location as 
the experimental curve. Here we see that the spectrum is beginning to resemble the experimentally determined one 
fairly closely, with its characteristic asymmetry. In the computed A± g spectrum from the single site operator, we find 
even more encouraging results. For a 16 site model, the SSO spectrum shown in Fig. || is peaked at almost exactly 
the same location (about 4.2J) as the experimentally-determined spectrum. The scaling used here is the same scaling 
used to match the peak heights for the B\ g spectrum. It is interesting to note that this scaling, chosen independently 
of the A\ g results, puts the peak at exactly the correct height. 

We find similarly encouraging results with a 26 site model. In Fig. [|, we again see the B\ g spectra for the SSO 
(solid curve of part (a)) and the Fleury-Loudon-Elliott operator (solid curve, part (b)) compared with experiment 
for La2Cu04 (dashed curves). Again we see that the SSO spectrum has more activity in better proportions than the 
Fleury-Loudon-Elliott curve does. The two-magnon peak (the maximum) is shifted slightly closer to what experiment 
shows, and there is more broadly distributed four-magnon activity in the SSO curve. In short, the asymmetry and 
line broadening seen in experiment is better suggested by the SSO spectrum than by the pure Fleury-Loudon-Elliott 
spectrum. If we again broaden the SSO spectrum and shift it slightly, in the same manner as for the 16 site model, we 
see in Fig. ^| that we have a better approximation of the experimentally-determined spectrum. In the A\ g spectrum 
shown in Fig. |^, we see the same encouraging signs of extra activity from a larger model, as compared with Fig. ^. 

Lastly, it should be pointed out that the goodness of the fit does appear to improve with increased system size. 
It would be helpful to see these calculations performed for larger systems. The SSO lacks translational symmetry, 
unfortunately, which prohibits many reductions in Hilbert space size that would otherwise be possible. For the 
moment, exact diagonalization of larger systems is beyond the capabilities of the authors' computing facilities. Only 
the Quantum Monte Carlo methodcJ can deal with substantially larger sizes and should prove specially informative. 

V. CONCLUSION 

In this paper, we have presented numerical data from exact diagonalization studies that suggest that improved 
understanding of magnetic Raman scattering in the insulating cuprates can result from a modification of the Fleury- 
Loudon-Elliott Raman operator. This assumes that phonons participate in the Raman scattering process, acting 
as momentum sinks and allowing for Raman scattering from two-magnon states with arbitrary total momentum. 
This provides a natural explanation for comparable Raman scattering in A\ g and B\ g configurations, and leads to 
a broadening of the spectra. This is achieved without invoking substantial modulations of local exchange constants, 
which can strongly effect long-wavelength properties. Due to limitations of sizes, the results presented are not fully 
conclusive about how close this brings the theoretical results to the experiments. Quantum Monte Carlo simulations 
may prove helpful in this regard. 

Given the large number of experiments on insulating cuprates, which can be modeled in terms of the square-lattice 
Heisenberg model, with a single nearest-neighbor exchange constant J, it seems natural that this be regarded as 
a good model for this system unless clear evidence to the contrary emerges. Raman scattering by itself cannot be 
invoked to justify more fancy terms such as ring-exchange terms for these systems. Raman scattering is also not 
the ideal ground for establishing the existence of spinons and other exotic excitations, although in the cuprates it 
definitely leaves room for such exotic physics. As more direct probes of quasiparticles, photoemission and neutron 
scattering can be more persuasive in this regard. 
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FIG. 1. (a) 16 site Fleury-Loudon-Elliott type Bi g scattering spectrum (solid) compared with experiment (dashed), (b) 
Single site operator scattering for the same system (solid) compared with experiment (dashed). 



FIG. 2. The SSO spectrum of Fig. |l|(b) broadened and shifted slightly (solid) to demonstrate the goodness of fit with 
experiment (dashed). 
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FIG. 3. The 16 site SSO A\ g spectrum compared with a sketch of the experimental data. Note that no shifting is necessary 
in this example. 

FIG. 4. (a) 26 site Fleury-Loudon-Elliott type Bi g scattering spectrum (solid) compared with experiment (dashed), (b) 
Single site operator scattering for the same system (solid) compared with experiment (dashed). 

FIG. 5. The SSO spectrum of Fig. ^(b) broadened and shifted slightly (solid) to demonstrate the goodness of fit with 
experiment (dashed). 



FIG. 6. The 26 site SSO A\ g spectrum (solid) compared with a sketch of the experimental data (dashed). 



7 







1 

(a) 

/ 

/ 

/ 

/ 

/ 

4 1 -J—h- 


1 i 

/ \ 
/ \ 
/ \ 
/ \ 

' \ 

\ 

~\ 

\ 

\ 

N 

L. 1 ^ 1 /A 1 h 


(b) 

/ 

/ 

/ 


1 \ 

1 \ 
1 \ 

\ 

A 

\ 

-\ 

V 



2 4 6 8 

co/J 



1 




1 



0.8 




